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ABSTRACT 



o 



In this paper we calculate the delay of the arrival times of visible photons on 
the focal plane of a telescope and its fluctuations as function of local atmospheric 
conditions (temperature, pressure, chemical composition, seeing values) and telescope 
diameter. The aim is to provide a model for delay and its fluctuations accurate to 
the picosecond level, as required by several very high time resolution astrophysical 
applications, such as comparison of radio and optical data on Giant Radio Bursts from 
optical pulsars, and Hanbury Brown Twiss Intensity Interferometry with Cerenkov 
light detectors. The results here presented have been calculated for the ESO telescopes 
in Chile (NTT, VLT, E-ELT), but the model can be easily applied to other sites 
and telescope diameters. Finally, we describe a theoretical mathematical model for 
calculating the Fried radius through the study of delay time fluctuations. 

Key words: atmospheric effects - methods: statistical. 
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1 INTRODUCTION 

Several high time resolution instruments are able to measure 
the arrival time of visible photons with an internal precision 
in the range 10 - 100 picoseconds, and refer those arrivals to 
the commonly used UTC scale with an accuracy of the o rder 
of 500 picoseconds ( see fo r instance, Barbieri et al. |2010| and 
Naletto & Barbieri. [ioogh . 

The step of referring the arrival times to UTC is usually 
done with the intermediary of time signals broadcasted over 
radio frequencies (e.g. GPS, GLONASS, Galileo GNSS), sig- 
nals which are very accurately corrected for atmospheric 
propagation effects like ionospheric scintillation and wet tro- 
posphere refraction. In the usual astronomical applications, 
the arrival times of optical photons are not corrected to the 
same accuracy. 

The possibility to perform such correction is actually shown 
by accurate laser ranging to geodetic satellites (recall that 1 
nanosecond in vac uum corresponds to 30 cm). For instance, 
Krai et al. (|2005l ) quote a precision of few picoseconds by 
taking into account the atmospheric seeing. Moreover a re- 
cent discussion has been performed by Dudy D.Wijaya and 
Fritz K. Brunner (2011) on the atmospheric range correc- 
tion for two-frequency Satellite Laser Ranging (SLR) mea- 
surements. 

Motivated by our own very precise time meas urements on 
celestial so urces with Aqueye and Iqueye (|20ld and Naletto 
fc Barbieri. |2009|) . we have undertaken the calculation of the 
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delay and delay dispersion of visible photon arrival times in 
the usual conditions prevailing in astronomical observato- 
ries. In the first step of our pr ocedure, the Marini-Murray 
model (Marini & Murray, 1 19731 ) is used to calculate the cor- 
rection Af? to the optical path of photons in air. 
Through this calculation, the atmospheric refractive index 
n and a fixed delay time independent of the photometric 
night quality is derived. Then the photon paths are corre- 
lated with the astronomical seeing. 

Through this relation we derive a statistical set of delay 
times as function of the Fried radius ro and telescope di- 
ameter. Finally, the difference of these delay times with the 
fixed delay gives the fluctuations. 

Reversing this model we also developed a theoretical mathe- 
matical model for the ro calculating through the observation 
of these fluctuations. 

The results are expounded in the present paper, which is 
organized as follows: 

• Section 2 describes the modified Marini-Murray model 

• In Section 3 the refraction index and delay time is cal- 
culated 

• In Section 4 the Fried radius is introduced and the see- 
ing effects on the image are recalled 

• In Section 5 the geometric and physical delay time fluc- 
tuations is introduced 

• In Section 6 we describe the model for the total delay 
time fluctuations 

• In Section 7 we calculate ro through the study of delay 
time fluctuations 
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2 MARINI-MURRAY MODEL FOR THE 
REFRACTION CORRECTION 

The Marini-Murray model is based on a n expa nsion of hy- 
pergeometric functions (Marini & Murray. 119731 ). The model 
relies on hydrostatic equilibrium and the barometric equa- 
tion. It considers a hydrostatic water vapor distribution and 
the water vapor behavior such as that of an ideal gas. The 
refraction correction is given as a function of temperature, 
pressure, vapor partial pressure, gravity acceleration and 
universal gas constant. The model also includes the values 
of water vapor molar mass and air average molar mass. Fi- 
nally the refraction correction is calculated in relation to 
the Earth radius, the site altitude and latitude. Through 
these parameters we can determine Optical Path Length 
(GPL) correction as described in detail below. According to 
the original Marini-Murray model, the refraction correction 
/S.R = Ri — R (see Fig. [T]), is calculated as: 



AR^ 



/(A) 



91+92+ 93 



91+92 + 93 



where gi, 32, 33 and g{()),H) are defined below, 6ui is 
angle of elevation, (j) is the latitude and H is the altitude of 
the observatory. The several terms are defined as: 
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In the original model the function /(A) was defined by 
with the following formula: 

r^^^r 0.0164 0.000228 
/(A) = 0.965 + + —^^ 

However, the refraction index formula used in the orig- 
inal Marini-Murray formalism is valid for a limited wave- 
length range (0.40/im 0.60/xm). Thi s is d ue to the model 
assumptions. For this reason, Ciddor (Il996l) devised a more 
refined model valid for a wider wavelength range (0.30/im-^ 
2.00/im). T his re fined model is described by the formula 
(Riepl et al. 120011 ): 



Ai? = fGr{\) 



91 +92 
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where the dispersion formula, taken from Ciddor (Il996l ) 
normalized to the wavelength of 0.6943/im, gives: 



fGr{\) = 
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Figure 1. Optical path length (OPL) R and its variation Ai? 
due to atmospheric refraction. R is the OPL in vacuum, R\ is 
the OPL in atmosphere, R^ is the Earth radius which the model 
calculates according to the latitude, 9zen is the Zenith angle and 
H is the site altitude. 



The values of the involved a-dimensional constants re- 
sulting from the normalized dispersion formula are: 
fco=238.0185, fci=205.0638, fc2=57.362 and ^3=5.944936. 
In the formulae, T [K] is the temperature, P [mh] is the to- 
tal pressure, and P^^ [mb] is the water vapor partial pressure, 
g — 9.784 [-^j is the acceleration of gravity at the equator. 
The gas constant is taken as Rg ~ 8314.36 ■ moZ] . The 
following molecular values are used: 



M- 

L mol J 

-2-1 

L mol J 



■ Water vapor 
Air average. 



• M„ = 18.016 

• Mi = 28.966 [; 

Equation [1] fills the Marini-Murray simplifications with 
a hybrid approach in which the dry and wet refraction delays 
are treated separately. 



3 REFRACTIVE INDEX AND DELAY TIME 
CALCULATION 

Through the refraction correction path AT? we can calculate, 
in a first approximation, the mean atmospheric refractive 
index n\. In fact, we note that the ni fiuctuation is of the 
order of 10~^ for 1 ps (see Section [5. 2|l . If we assume that 
ni has a constant distribution: 



m = 



O PL + AR 
OPL 



(3) 



where the geometric optical path length (OPL) is given 
by the formula: 
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OPL = 



AA-SA 



COS(ti^en 

where: 



(4) 



• AA=Atmosphere indicative thickness (In this paper we 
have considered the AA = lOKm) 

• SA=Site altitude 

The refractive angle Oi is then calculated as: 

6i = arcsin ( —\ 
\niJ 

With the value of the refractive index ni the delay time 
ti is calculated through the formula: 

OPL OPL 



ti 



ni 



c/m 



(5) 



where c is the 
(c = 299792.458 [^]). 



velocity of light in vacuum 
Equation [5] allows us to cal- 
culate the propagation delay time in the atmosphere as a 
function of the refractive index n, including the geometric 
effects. 

The first term in Equation [5] is the ratio between the 
geometric OPL and the speed of light in the atmosphere. 
The main contribution of the delay time comes from the 
denominator of this ratio. 



3.1 Application of the model to ESO astronomical 
sites for the delay time 

In this Section we calculate the delay time for the three 
Chilean sites of ESO telescopes. Table [T] shows the sites 
characteristics. We assume an average ground temperature 
of 288 K and an average ground relative humidity of 20% 
for the three sites. Table [5] shows the obtained values from 
the simulation (through the Equation [5]) at different Zenith 
angles. The delay time is calculated as a function of wave- 
length (A — 0.632^m) and has been calculated for a zenith 
angle ranging from 0° to 60° . The data show the increase of 
the delay time with the 9zen following the trend of the cos 
function. We note that this delay does not depend on the 
telescope diameter and has variations in the order of 10~'^ns 
between Paranal and La Silla. These variations are due to 
the difference in altitude and geographical coordinates of the 
sites. 



4 CALCULATION OF THE FRIED 
PARAMETER AND SEEING 

The atmospheric optical turbulence introduces variations on 
ti according to refraction fluctuations. We now recall one 
of the main parameters of the astronomical seeing, namely 
the Fried's radius ro, which defines the average size of the 
turbulent cell. This parameter will be useful to introduce the 
concep t of p hotons delay time fluctuations. Fried has shown 
f Fried. Il965l ). within the limits of validity of the Kolmogoroff 
law, that ro is expressed by the formula: 



ro = 



0.423 



4TT-' 







C„ ■ dz 



(6) 



Table 1. Geographic Characteristics of the Sites. 



Site 


LAT. 


LONG. Altitude 


Telescope Diameter 








Km 


m 


La Silla 


-29°15 


-70° 43' 


2.347 


3.58 


Paranal 


-24° 37 


-70° 24' 


2.630 


8.20 


Armazoncs 


-24° 35 


-70°11' 


3.064 


42.00 


Table 2. Delay time 


vs Zenith angle 


for different sites (A = 


0.632^m). We note that this delay time does not depend on 


the telescope diameter and varies slightly with the location of 


the site. 










Zenith Angle 


n 


15 


30 


45 60 


La Silla 


(ns) 


15.53 16.24 


18.58 


23.14 30.72 


Paranal 


(ns) 


15.52 16.23 


18.57 


23.13 30.71 


Armazones 


(ns) 


15.52 16.23 


18.57 


23.13 30.71 



80 ■ 10^ 



r 



and the temperature structure parameter C|.(a;) is de- 
fined through the formula: 

2 _ [T{x)-T{x + Ax)f 

— 2 

this parameter is expressed in (°C)^-m~^ and expresses 
the temperature variations between two locations at a dis- 
tance Aa;. 



4.1 Seeing effects on the images 

The seeing produces scintillati on, sm earing and motion of 
the image. Roddier (Roddier, Il98lh has obtained the fol- 
lowing approximate expressions for the calculation of these 
three effects. 



4-1.1 Scintillation 

The image scintillation, as a function of Cf^, in approxima- 
tion is given by the following formula: 



1 



/ (Cos(e.en))^ 



Cn (z) ■ ■ dz 



where is the refractive index structure parameter: 



where D is the diameter of the telescope. Scintillation 
however does not affect the dispersion of the arrival times, 
and therefore is not taken further into account. 



4-1.2 Image Smearing 

The light from the point source is spread over an area having 
a Full Width Half Maximum (FWIfM) given by: 

FWHM = 0.98 — 
ro 

The value is in arcsec- The amplitude of this effect is 
independent of the pupil diameter. 
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4-1.3 Image Motion 



The motion of the image, as a function of A, telescope di- 
ameter (D) and ro, is given by: 



Source 



u'^{x) = a^{y) = 0.18 • • D~3 ■ r~ 



(7) 



This motion is in arcsec, and can be expressed in linear 
units (meters) by means of the formulae: 



Ax = ni • OPL ■ sin^ a'^(x) 

Ay — ni ■ OPL ■ sin^J cr^{y) 

The latter effect is very important for the fluctuation 
of delay times. One can calculate a new optical path as a 
function of this image motion, and then a new delay time 
(See Figure [5]). Figure [5] shows the optical path length in 
vacuum and the graphical representation (not to scale) of 
its variations. 

We define OPLi = ni ■ OPL and OPL2 = n2 ■ OPL, where 
ni is calculated using the Marini-Murray model (see For- 
mula [S} and OPL through the Formula |4l The difference 
between the previously calculated time and this new one 
gives the fluctuation. The total motion is given by the for- 
mula: 



y 



(8) 



This motion is << OPL, then we can calculate the new 
OPL through the formula: 



(9) 



OPL2 = ^/OPLl + A2 

where: 
OPLx = m ■ OPL. 



5 GEOMETRIC AND PHYSICAL DELAY 
TIME FLUCTUATION 

In this Section we calculate the delay time fluctuation. It is 
due to two contributions: one correction is due to the change 
in photon path (geometric), and the other (physical) is due 
to the fact that the photon is traveling in a medium where 
the refractive index is changed. The delay time (t2,G(»'o, D)) 
due only to the geometric variation of the OPL is given by: 



ni ■ OPL2 OPL 



Through this new time we get the geometric fluctuation 
of the delay time Atciro, D): 



AtG{ro,D) = \t2,Giro,D)~ti\ 
Taking into account that: 

m ■ OPLi OPL 
ti = 

c c 

we obtain the following formula: 

Atciro, D) = — ■ \OPL2 - OPLi I . 
c 



(10) 



(11) 
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Figure 2. Schematic diagram of the OPL variation due to the 
image motion. Figure shows the optical path length in vacuum 
and the graphical representation (not to scale) of its variations. 
We define OPLi = m ■ OPL and OPL2 = n2 ■ OPL. Considering 
the motion of the image cr < < ro we can assume that relevant x 
and y refractive index gradients are unexpected. 



5.1 Physical delay time fluctuation 

In this Section, we calculate the delay time (t2,p('"o, L))) due 
to the physical variation of the atmosphere. In the previous 
section we saw how the OPL changes due to the atmosphere 
(image motion). The change in OPL also induces a change 
in the refractive index. In fact, fluctuations in optical path 
length must correspond to refractive angle variations. We 
can calculate the refractive angle fluctuation assuming the 
atmosphere uniform distribution: 



A6I : 



OPLi 



^OPL\ + A2 



(12) 



The corresponding atmospheric refractive index n2 , ob- 
tained by Snell's law, is: 

n\sinQ\ 



n2{ro,D) 



sinidi + A6) 



(13) 



This refractive index is a function of ro and the tele- 
scope diameter(L>) . 



5.2 Calculation of physical delay time fluctuation 

With the value of the refractive index n2 (Equation I13[) we 
can calculate a new physical delay time as a function of the 
ro and the telescope diameter (D), through the formula: 



t2,p{ro,D) = 



n2 ■ OPL _ OPL 

c/n2 c 
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Fluctuation vs. ro 




Figure 3. Fluctuation delay time as a function of the ro- This 
simulation pertains to La Silla (A = 0.632/im, Zenith angle = 
10°). 



Fluctuation vs. Zenith Angle 

(Diameter Variation) 



40 - 



D=3.58 m 

D=8.20 m 
— D=42.00 m 




10 20 30 40 60 60 70 

Zenith Angle (°) 

Figure 4. Delay time fluctuation as a function of the telescope 
diameter and the Zenith angle variation (A = 0.632/im, ro = 
15cm). The simulation compares La Silla (NTT), Paranal (VLT) 
and Armazones (ELT). 



Table 3. Fluctuation vs Zenith angle for different sites (A : 
0. 632^771, ro = 15cm). 



Zenith Angle 


(°) 





15 


30 


45 


60 


La Silla (3.58m) 


(ps) 


9.8 


10.4 


12.4 


16.9 


27.5 


Paranal (8.20m) 


(ps) 


8.2 


8.7 


10.4 


14.2 


23.1 


Armazones (42.00m) 


(ps) 


5.8 


6.1 


7.3 


9.9 


16.2 



where OPL is given by the Formula [l] Then, as in 
Section [S] we obtain the fluctuation due to this variation 
Atp{ro,D): 



Atp{ro,D) = |t2,p(ro,D) 
where: 

ni ■ OPLi OPL 



(14) 



we obtain: 



^ , OPL I 2 
Atp{ro,D) = \nl 



2 I 



(15) 



Through the Formula [15] we can estimate that a delay 
time fluctuation of 10.0 ps corresponds to a refractive index 
variation An of 3.3 x 10~*. This gives an idea of the error 
propagation. 



6 CALCULATION OF DELAY TIME 
FLUCTUATIONS 

Consequently, the total fluctuation of the delay time is the 
sum of the geometric and physical component: 

Atp = Ata + Atp 

In this sum the physical component is larger than the 
geometric component of four orders of magnitude. Finally, 
we express this fluctuation through the extended formula: 



AtF{ro,D) 



OPL 



nl ■ OPL 



msinOi 



^/2s 



+ 



sin{ei + A9) 



(16) 



where: 



AO : 



2sin \ K ■ ■ D 



Substituting the constants with these values (Roddier, 

Il98ll ) K = 0.18, a = 1, /3 = 3, 7 = 5, 77 = 3 and e = 2 we 
obtain the final formula used in the model: 



AtF{ro,D) = 



+ 



OPL 

c 

where: 



A9 = arcos 



nl ■ OPL 



nisinOi 



V2sin\ 0.18 ■ A2 ■ D-3 . r," 



sin{9i + A9) 



2 



(17) 



We note that this delay time fluctuation depends on the 
wavelength and on the telescope diameter. In particular, it 
also depends on the Fried radius. Figure [3] shows the trend 
of this fluctuation as a ro function. The simulation is done 
for La Silla (A = 0.632/im, Zenith angle = 10°). 

6.1 Application of the model to ESO astronomical 
sites for the delay time fluctuations 

In this Section we apply the previously described model to 
the three Chilean sites of ESO telescopes. We present the 
results obtained through the Formula 1 171 Table [3] shows the 
simulation results for each site, in particular the fluctuations 
variation as a function of the telescope diameter. We note 
that the fluctuation decreases with a larger diameter. In this 
case, between La Silla and Armazones telescopes the delay 
time fluctuations decreased of about 40%. Figure [4] shows 
the results of this model for each site as a function of Zenith 
angle. 

Table |4] shows the simulation results for three wavelengths. 
In this case the simulation is done for La Silla. Figure 
[5] shows the Table [4] results and Figure |6] is the zoom 
of these trends. In fact, we note that in the range A = 
0.416 — 0.799/im the fluctuation changes only by 4%. 
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Fluctuation vs. Zenith Angie 

(Wavelength Variation) 




30 40 

Zenith Angle (°) 

Figure 5. Delay time fluctuation as a function of the wavelengtli 
and tile Zenitli angie variation. Tlie simuiation pertains to La 
Siiia {D = 3.58m, ro = 15cm). 



Table 4. Fluctuation vs Zenith angle for different wave- 
lengths, La Silla {D = 3.58m, tq = 15cm). 



Zenith Angle 


n 





15 


30 


45 


60 


A = 0A16^lm 




10.1 


10.7 


12.7 


17.4 


28.3 


A = 0.632^m 




9.8 


10.4 


12.4 


16.9 


27.5 


A = 0.^99^lm 


(ps) 


9.7 


10.3 


12.2 


16.7 


27.3 



7 INVERTING THE MODEL TO CALCULATE 
FRIED RADIUS 

In this Section, we reverse the Formula ll7l in order to derive a 
mathematical expression giving the ro value from the study 
of the delay time fluctuations (see Appendix The ro 
value is expressed as a function of the telescope diameter 
(D), of the wavelength (A) and of the variable (x)- This 
variable is obtained through the observation of the delay 
time fluctuations and the ni calculation. The model gives 
us the following formula for the Fried radius calculation: 



ro = 0.54 ■ 



A 5 • arctg ^ (xi.a) 



(18) 



where: 



Xl,2 



Z - tg^Oi 



Fluctuation vs. Zenith Angle (Zoom) 

(Wavelength Variation) 



---A=0.416 

1=0.632 
-^1=0.799 



O 

is 



10,0 




14 16 

Zenith Angie ("[ 



Table 5. Simulation of tlie ro calculation through the delay 
time fluctuations for La Silla (A = 0.632/im, telescope diame- 
ter D = 3.58m). 



Delay time fluctuation Fried radius 


Atp{ro,D) 


ro 


ps 


cm 


27.0 


5.0 


18.2 


10.0 


10.2 


20.0 



Figure 6. Zoom of FigurefS] 



and: 

c- AtFiro,D) +nl ■ OPL 
n\ ■ OPL 

Through this model we have a new way to calculate the 
Fried radius as a function of the OPL and average refraction 
index of the atmosphere. Table [5] shows a simulation for La 
Silla, the mathematical details are given in Appendix [Al 



8 CONCLUSION 

We have shown in this paper that the delay and dispersion 
introduced by the terrestrial atmosphere in photon arrival 
times is significant to the level of tens of nanoseconds and 
tens of picoseconds respectively. 

These values mean a severe degradation of the time tagging 
performances of modern astronomical detec tors such as the 
Avalanche Photodiode (e.g. Barbieri |201C| and Naletto & 
BarbieriJOOO). Furthermore, they imply that the very accu- 
rate timing signals available e.g. via radio signals by GPS or 
other Satellite Navigation System are not exploited to their 
full extent in the astrophysical field, as instead is done in 
geodetic satellite ranging applications. The scientific results 
of very high time resolution astrophysics can therefore be 
degraded by such neglect. In particular, the E-ELT will pro- 
vide a 25-fold increase of photon fiux over existing telescopes 
thus opening the way to 'quantum' astronomy (e.g. Dravins 
et al.^006). The algorithms expounded in the present paper 
overcome these limitations. 

In Section [7] we have described a theoretical mathematical 
model for calculating the Fried radius through the study of 
delay time fiuctuations. Table[5]shows the simulation results 
for La Silla. This is a completely new method for the study 
of atmospheric turbulence, in future work we want to corre- 
late the delay time fiuctuations experimental data with the 
traditional techniques for the Fried radius calculation. This 
will allow us a broader view of Earth's atmosphere and its 
infiuence on propagation of photons. 
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comes negligible. 

Then we consider the formula: 



AtF{ro,D) 



OPL 



riisinBi 



i{9i + A9) 



2 



(A2) 



we isolate the trigonometric functions: 
AtF{ro,D) +nlOPL _ 



sm Oi 



n\OPL 
we substitute: 



sm2(6ii + AS) 



c- MF(ro,D) +nl ■ OPL 
nl ■ OPL 



(A3) 



Using the trigonometric addition formulas and param- 
eters of the sin and cos functions we get: 

[Z - tg^ei)x^ + 2- Z ■ tgOix + Z ■ tg'^di - tg'^9i = 

where: 

X = tg^e (A4) 

and solving in the same variable we obtain: 

-tg9i ± y/tg*6i{Z- l)+tg'^ei ■ Z 
-^''^ " Z - tg^9i 

We now consider the equations lAll and IA4I obtaining 
the following relationship: 



A6 = arctgx = arcos 



2sin V0.18A2L)-3r„ 



(A5) 



Finally, from equation |A5] we explicitly the ro value and 
we obtain the following formula: 



ro = 



0.54 ■ D 5.^5. arctg s (xi,: 



(A6) 



APPENDIX A: INVERTING THE MODEL 



We start from Equation 1171 relating the Atpljo, D) to the 
Fried parameter ro: 



AtF{ro,D) 



nl ■ OPL 



+ - 



OPL 



n\sin6i 
sin(ei + AO) 



V2sin\/ 0.18 ■ A2 • D-t . 2 - 1 



2 

ni 



+ 



where: 
A6 = arcos 



2siny0.18 ■ A2 ■ D'a • + 1 



(Al) 



In the case of studying the delay time fluctuations for 
the ro calculation, the geometric delay time fluctuation be- 
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